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Abstract 

Recent theoretical approaches to community structure and dynamics reveal that 
many large-scale features of community structure (such as species-rank distributions 
and species-area relations) can be explained by a so-called neutral model. Using 
this approach, species are taken to be equivalent and trophic relations are not taken 
into account explicitly. Here we provide a general analytic solution to the local 
community model of Hubbell's neutral theory of biodiversity by recasting it as an urn 
model i.e. a Markovian description of states and their transitions. Both stationary 
and time-dependent distributions are analysed. The stationary distribution — also 
called the zero-sum multinomial — is given in closed form. An approximate form 
for the time-dependence is obtained by using an expansion of the master equation. 
The temporal evolution of the approximate distribution is shown to be a good 
representation for the true temporal evolution for a large range of parameter values. 
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1 Introduction 



Understanding the global patterns of biodiversity and its dynamics at different time scales 
remains a great challenge for ecological science (Rosenzweig, f995; Wilson, 2003). One 
of the key features that defines community structure is the relation between range and 
abundance. How the community structure develops in time and how species are spatially 
distributed largely define the field of macroecology (Brown, 1995). In this context, an 
important step to unify biogeography and biodiversity has been achieved by Hubbell 
(Hubbell, 2001; Bell, 2001) through the formulation of a neutral theory. 

The mathematical framework employed by Hubbell allows for speciation processes to 
be integrated with the MacArthur- Wilson theory of island biogeography. In this way, the 
neutral theory predicts some universal features that can be tested by direct analysis of 
species-abundance distributions and other large-scale measurements. In Hubbell's theory, 
two key quantities largely determine the steady-state distributions of species richness 
(as well as relative species abundances on local and large geographic scales). These 
two parameters are the so-called biodiversity number and the immigration (dispersal) 
rate. Under Hubbell's assumptions, the ecological properties of every individual in the 
population are assumed to be identical. 

In a neutral model of this type, individuals compete for the same pool of resources, 
but chance events are responsible for the identity of the final winner (s). The dynamics of 
each species is thus path-dependent and a Markovian description of their time evolution 
is appropriate. Under the assumption of a balance between birth, death and immigration 
rates, the neutral theory is able to reproduce the quantitative patterns of species distri- 
butions that are well known from the ecological literature. It also permits the generation 
of several nontrivial and testable quantitative predictions about biodiversity and biogeog- 
raphy. In particular, the theory predicts that rare species would typically be recent in 
terms of their origination. In relation to conservation biology, a neutral community in 
which species are essentially equal would be very fluid, with frequent replacements. If 
true, protected areas should be larger than those expected for stable communities with 
species closely adapted to given niches. 

Formally, Hubbell's theory is the ecological analog to the neutral theory of genetic drift 
in genetics (Kimura, 1983; Ewens, 1972; Karlin and McGregor, 1972). Early attempts 
to incorporate the neutral approach from population genetics (Caswell, 1976; Hubbell, 
1979) mainly highlighted the relevance of drift in community dynamics, providing evi- 
dence for a global view of ecosystems in which competitive forces, ecological niches, and 
even trophic interactions could be ignored in the pursuit of a better understanding of 
biodiversity dynamics. More recent work incorporated these ideas in an explicit way 
(Hubbell, 1997; Sole and Alonso, 1998) and Hubbell's recent book provides an extensive, 
unifying account of these (Hubbell, 2001). The starting point of neutral models is a ran- 
dom community that evolves towards an equilibrium state under a stochastic birth-and 
death process incorporating dispersal. At high immigration rates, Hubbell's theory pre- 
dicts a logseries distribution for the abundance of species in the local community, while 
when the immigration coupling between the metacommunity and the local community is 
lower, a lognormal-like shape is obtained for this distribution. Within Hubbell's approx- 
imation, these distributions are shown to be particular cases of what he denotes as the 
zero-sum multinomial (Hubbell, 2001). 

Hubbell's model for local communities is similar to that proposed in Sole et al, (2000) 
and analysed in McKane et al, 2000. There we took advantage of a mean field argument 
to find an analytical form for the stationary distribution for the probability of finding 
species having an abundance of n individuals. In addition, we studied in detail its time 
behaviour using different approximations. Furthermore, our simplified approach based 
on this mean field argument allowed us to recover the scaling relationship between the 



2 



fraction of links actualised and the number of species in a community — the so-called 
C*-S relation, and gave conditions in which such a relation arose. 

Within Hubbell's mathematical framework the dynamical stochastic models were nu- 
merically solved and the equilibrium properties analysed. In this paper we present an 
analytic, general solution of Hubbell's model for the local community dynamics, that 
provides the stationary species-abundance distributions together with the time evolution 
from the initial state towards the stationary distribution. 



2 Formulation of the theory 

Hubbell's theory concerns populations on two scales: local communities and regional 
metacommunities. To explain the model and derive the equations in the simplest possible 
way, we will use the language of urn models (Feller, 1968; Johnson and Kotz, 1977). This 
is a natural description when the stochastic dynamics in one time step only depends on the 
state of the system at the beginning of the time step (in other words is a Markov process). 
It also provides us with a concrete picture of the process which aids the derivation of the 
governing equation for the model. 

We begin by considering the model in a limit where the two levels of description are 
uncoupled. This allows us to focus only on the local community. We assume that there 
are N individuals of species % in the local community, with the total number of individuals 
of all species being J, that is, J = Yh=i Ni where r is the total number of species. The 
model is defined by picking one individual at random from the local community, killing 
it, and then replacing it by an individual also drawn from the local community. In terms 
of the associated model this corresponds to having N balls of colour i (i — 1, . . . , r) in 
the urn. If we focus on one particular colour, j, the probability that the number of balls 
will decrease from Nj to Nj — 1 during one time step is 



since a ball of colour j must be discarded and one of any other colour replaced for such 
a transition to occur. On the other hand, the probability that the number of balls will 
increase from Nj to Nj + 1 requires that a ball of any other colour but j must be discarded, 
and one of colour j be replaced. Therefore 

^(^, + 11^) = ^^-^. (2) 

The whole point of the model, however, is to couple local communities and regional 
metacommunities. This is achieved by choosing a replacement ball from the urn only (1 — 
m) of the time. For the rest of the time it is chosen from outside the urn. The probability 
of picking a ball of colour j from this external source is defined to be Pj, and corresponds 
to assuming that the replacement individual comes from the regional metacommunity 
where species i has a relative abundance of Pj. The transition probabilities (0) and (|2*| 
now read 

W(N j -l\N j ) = {l-m)^y0& +m^(l-P j ) (3) 

W (Nj + l\Nj) = (1 - m) {J ~ Ni) + m ^'f^ Pj . (4) 

The change in the probability that there are Nj balls in the urn from time t to the time 
after one time step has elapsed consists of four contributions. Two of these correspond 



and 
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to an increase in this probability (due to transitions from (Nj — 1) and (Nj + 1) to Nj) 
and two to a decrease (due to transitions from Nj to (Nj + 1) and (Nj — 1)). The balance 
equation showing this change is: 

AP(Nj,t) = W(Nj\Nj-l)P(Nj-l,t) + W(Nj\Nj + l)P(Nj + l,t) 

- {W(Nj + l\Nj) + W(Nj-l\Nj)}P(Nj,t). (5) 

Compared with the long time scales we are interested in — during which many transi- 
tions will take place — the step size is very small, and we may take the limit in which 
AP(Nj,t) — > dP(Nj,t)/dt. The resulting equation is a master equation for the probabil- 
ity P(Nj, t) (Van Kampen, 1981; Gardiner, 1985). Some care is needed with the boundary 
conditions on this equation: clearly the cases Nj = and Nj = J are special cases since 
there can be no transitions which reduce Nj in the former case or which increase Nj in 
the latter case. One possibility is to write two separate equations for these special cases. 
However there is no need for this if we first observe that some of these conditions are nat- 
ural consequences of the form of the transition probabilities. For example, the expressions 
in and (J3J) are both zero if Nj = and Nj = J respectively. So as long as we agree 
to impose the formal definitions W(0\ — 1) = and W(J\J + 1) = the same master 
equation may be used for all states. In addition, an initial condition needs to be imposed 
to complete the specification of the problem. Typically, the number of individuals in the 
local community at t = will be given: P(Nj, 0) = 5jv-,jv - 

The mathematical formulation of Hubbell's theory described above can be directly 
mapped on to another dynamical model of a multispecies community which we introduced 
a few years ago (Sole et at, 2000; McKane et al, 2000; Sole et al, 2002). In this case 
though, the nature of the interaction depends on the "score" between one species and 
another, and a form of mean field theory had to be used in order to describe the dynamics 
by such a straightforward dynamics. In terms of the notation we have used above — iV 
denoting the number of individuals of a particular species and J denoting the total number 
of individuals of all species — the transition probabilities of this model are (Sole et al, 
2000; McKane et al, 2000): 

and 

W(N -1\N)= C*(l - rfj + |(5 - l)j . (7) 

Here \i is the fraction of the time that replacing of one species by another can happen 
by chance, and not because the replacement individual belongs to a species which has a 
positive score against the first. It clearly maps into m. The other constants are S, the 
number of species, and C*, a parameter related to the degree of connectivity of the matrix 
of scores between the species. The precise form of the mapping is C* = 1 and Pj = S~ l . 

Since we have analysed this model extensively (McKane et al, 2000) we may simply 
deduce expressions for quantities of interest in the Hubbell theory by setting C* = 1, S = 
Pf 1 and jj, = m. 



3 Stationary state 



The most straightforward questions we can investigate concern the nature of the stationary 
state of the theory. Let us begin by introducing the abbreviations 



W(Nj 



s (J - Nj) 
1-m)^- 3 — +m(l 

J 1 



Pi 
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and 



g N =W{Nj + l\N. 



3) 



J 



. N 

m) h mPj 

J —1 3 



(9) 



The master equation now reads 

dP ^ l) = r Nj+ i P(Nj +l,t) + g Nj -! P(N, - 1, t) - {r Nj + g Nj } P(Nj, t) . (10) 

The stationary probability distribution, P s (Nj), is determined by setting dP(Nj)/dt = 0. 
This gives 

r Nj+1 P S (N, + 1) - g Nj P s (Nj) = r N . P a (Nj) - g N ,-i P,(A^- - 1) . (11) 

This is true for all Nj, which implies that P s (Nj) — gNj-i Ps(Nj — 1) — I, where / is a 
constant. Applying the boundary condition at Nj = 0, we find that 1 = and therefore 

r Nj+1 P s (Nj + 1) = g Nj P s {Nj) ; Nj = 0, 1, J. (12) 

To solve this equation, let us first assume that m^O. Then the rjv,- and g^, given by (jHJ) 
and Q are all non-zero and we can solve (|T2"j) by iteration to obtain 

r N . r N] _ l . . . n 

The constant P s (0) can be determined from the normalisation condition 

E p s(Nj) = P s (0) + E P.{Nj) = 1 • (14) 

Nj=0 Nj>0 

To simplify the algebra let us introduce some new notation for various combinations of 
parameters which naturally appear in the solution of the model. We write the transition 
probabilities as 

(1 — in) 



J(J- 



~k Nj {Nj — Nj) , (15) 



and 



where 



p . = P and N , = ?J^\ _ P , (17) 

J (1 — m) J V 1 - mj J 

Substituting the expressions (fT5|l and (JTHj) into (|T^|) gives an explicit representation for the 
P s (Nj) in terms of P s (0). An expression for P s {0) itself can be obtained by performing the 
finite sum which appears in ()14|1 . This sum can be performed analytically using properties 
of Jacobi polynomials (Abramowitz and Stegun, 1965). Alternatively, the mapping into 
the model defined by © and (J7J) can be used since the result for the P s is known in this 
case (Sole et al, 2000; McKane et al, 2000). One finds (see McKane et al, 2000, for 
details of the derivation): 



s{ j) \Nj) (3(P*,N*-J) 
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where (3(a, b) = T(a)T(b)/T(a + b) is the beta-function. 

It is interesting to note that in the case m = 0, where the local community is decoupled 
from the regional metacommunity, go = 0, and so from (|12|). since r\ 7^ 0, it follows that 
P s (l) = 0. In fact, since r Nj ^ for < Nj < J, we see from (JEJ) that P s (Nj) = for all 
< Nj < J. So with no interaction with the regional metacommunity, species j either 
disappears or becomes the only species there is in the local community. Therefore some 
degree of coupling is vital for biodiversity. 

In Fig 1, we have computed the stationary distribution for different parameter values 
and sizes of the system. The relative species abundance distribution predicted to occur 
in local communities — the zero-sum multinomial — by the unified theory of Hubbell can 
be readily computed even for high community sizes using the analytic formula (j!8|) . 



4 Time dependence 

Together with universal features displayed by the stationary patterns observed in mature 
communities, some common features are also observed when looking at how diversity de- 
velops in time. When an empty field starts to be colonized by immigrant species a new 
community gets formed and a pattern of species replacement develops. The transition 
from abandoned field to mature forest is one of the best known examples of ecological 
sucession and is common in many places after the abandonment of agricultural land. In 
temperate climates, a mature forest is the end point of sucession, in spite of the difer- 
ent potential initial conditions. The path towards the steady species-ranks distribution 
seems to be common to many different ecosystems (Hubbell, 2001). Furthermore, natural 
systems are continuously perturbed; any disturbance resumes the process of ecological 
succession. It is thus natural to ask: what predictions about this process can be made in 
the context of Hubbell's neutral theory? 

In the last section it was shown that a closed form expression could be obtained for 
the probability of finding Nj individuals of species j in the local community when the 
systems has reached the stationary state. In addition to this, just mentioned, we also 
wish to know how the community is assembled from a given starting point. This requires 
us to solve for the time-dependence of the model. It is not possible, in general, to carry 
this out exactly, since the transition probabilities (j!5|) and (JTfij) are nonlinear functions of 
Nj. It is nevertheless possible to get a very good approximation to P(Nj,t) by using the 
fact that in cases of interest J will be large. The approach which we will use, due to Van 
Kampen (1981), is rather technical and has been discussed elsewhere in some detail (Van 
Kampen, 1981; McKane et al, 2000), but the basic idea is quite simple. Therefore, we 
will avoid these complications, and quote relevant results using the correspondence with 
the transition probabilities © and (J7J). 

The key idea is to expand about the deterministic version of the theory. In the limit 
where the number of individuals becomes infinite, all stochasticity is lost, and the system is 
completely described by a deterministic equation. This equation is not known a priori, but 
if it can be established, an expansion in powers of J~ l could perhaps be set up to calculate 
corrections to the deterministic result which would be valid for large, but finite, J . Quite 
generally we would expect a plot of P(Nj,t) against Nj for fixed t to be approximately 
Gaussian for large J. The motion of the peak of this distribution would move with 
t according to the deterministic equation. Van Kampen's large J expansion gives the 
deterministic equation as the zeroth order ( J — > 00) result, with the next to leading order 
result giving a Gaussian distribution peaked at this value. Higher order contributions give 
corrections to this distribution, but they are usually so small for large J that they are of 
very little interest. Since a Gaussian centred on a given value is completely determined 
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by its width, there are only two things to find: (i) the deterministic equation, (ii) the 
width of the distribution. 

In practice one writes Nj = J<f>j{t) + J x l 2 x^ where <j>j(t) = \imj^ oc (Nj / J) is the 
fraction of j species which are present in the local community at time t in the deterministic 
limit. The variable 

Xj = -J= {Nj - J0,(t)) 

characterises the fluctuations away from the deterministic theory. We require 4>j(t) and 
(x^) ((xj) = 0). Using the correspondence between the two models we obtain (McKane 
et al, 2000) 

^ = m (Pj - , (19) 

where r = tj J is a rescaled time. This equation is easily understood: if (f>j is less than the 
abundance of species j in the regional metacommunity, then it increases. If it is more, 
then it decreases. The equation is easily solved to give 
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{r) = cj>0)e-™ + P j (l-e™ ). (20) 

Initially we ask that Xj(0) = 0, which means that 4>j(0) = Nj(0)/J = Nj^/J. Going back 
to the t variable gives 

<^(t) = ^2 e - m '/ J + Pj (1 - e - mtlJ ) . (21) 

In Hubbell (2001, Chapter 4), an alternative discrete-time formulation of this local 
community model is given. Obviously, both time discrete and time continuous formula- 
tions give rise to the same equations for the deterministic model counterpart (Hubbell, 
2001, page 110). However, he does not address the stochastic time-continuous formula- 
tion. Here we show that insight can be gained by finding approximate solutions to the 
time-dependent model. 

The width of the distribution is given by 



e~ 2mr 



) — 777 

+ A~—^{\ -2P j )e~ mT 
J m 



e -mr 



- 2(1 - m)A]re- lmr , (22) 



where Aj = (Nj t0 / J) — Pj. We have already commented that the probability distribution 
is a Gaussian to the order we have been working. Specifically, in terms of the quantities 
calculated above, 

where <f>j(t) and (xj) T are given by equations (J2~T|) and (j2*2|) respectively. 

In Fig. 2, we show the temporal evolution for P(Nj, t) computed both using a Gaussian 
approximation (Eq. (|23p) and the numerical integration of the master equation. The good 
agreement which is obtained is a reflection of the fact that community sizes J are taken 
to be large enough so that further terms in the large J-expansion are negligible. However, 
if the final stationary distribution does not have a Gaussian shape, more terms should 
be included in the expansion so as to capture the true temporal behaviour of P(Nj,t). 
Notice that, while the approximation given by Eq. (|2*3*|) is always represented as dotted 
or punctuated curves, in some cases these are not visible because they match the exact 
distribution so completely. 
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5 Conclusion 



The main aim of this paper has been to show that aspects of Hubbell's neutral model 
of local community biodiversity dynamics can be solved for exactly, and even if this is 
not possible, calculational schemes are available which provide very good approximations 
to the solution. Specifically, we have shown that the stationary properties of the model, 
which can be obtained from the zero-sum multinomial, can all be found exactly So, for 
instance, the mean value and variance of the number of individuals of species j, can be 
obtained from this probability distribution. The nature of the time evolution cannot be 
determined in closed form, but a controlled approximation based on assuming that the 
total number of individuals of all species, J, is large, is possible. This is an excellent 
approximation in most cases of interest, and we would expect that the results that we 
have obtained will be relevant in these situations. The applicability of our approximation 
scheme was checked by carrying out the numerical integration of the master equation (Eq. 
HU|) . The results, displayed in Fig. 2, confirm our expectations. 

While the results which we have reported describe the essential aspects of the solution 
of Hubbell's model, there are many other interesting features which are also amenable 
to analysis and for which definite, and well-controlled, results may be obtained. The 
structure of the metacommunity and the form of the colonisation curve are examples. 
These, and related questions, are presently under study, and we hope to report our results 
in a future publication. 
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Figure captions 



1. Zero-sum multinomial distribution. The analytic formula ()18|) has been used to 
compute the stationary distribution, P s (Nj), for different values of the abundance 
of species j in the metacommunity, the total number of individuals J and the prob- 
ability of immigration from the metacommunity, m. We have dropped the subscript 
j, which labels a particular species, in the figure. 

2. Temporal evolution of the probability, P(Nj,t), of having the j-th species repre- 
sented by Nj individuals. The temporal evolution has been computed using both 
the Gaussian approximation and the straightforward numerical integration of the 
exact master equation. In both cases, the initial number of individuals of the focus 
species was 0.8 x J. The relative abundance of the focus species in the metacom- 
munity was Pj = 0.1 also in both cases. We have dropped the subscript j, which 
labels a particular species, in the figure. 
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